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Abstract 


A discussion is given of ongoing work related to development and utilization of a numerical model for 
treating the fluid dynamics of ink Jets. The model embodies the complete nonlinear, time dependent » 
axi-symmetric equations in finite difference form. An earlier work treated continuous Jets in which periodic 
boundary conditions allowed study of local capillary Instability to drop formation in a moving reference 
frame. The present study includes the Jet nozzle geometry with no-slip boundary conditions and the existence 
of a contact circle. The contact circle is allowed seme freedom of movement, but wetting of exterior surfaces 
has not yet been addressed. The principal objective in current numerical experiments is to determine what 
pressure history, in conjunction with surface forces, will lead to clean drop formation. 

I. Introduction 


Recent experimental fluid drop studies related to ink Jet printer design have provided an extensive input 
to guide development of a comprehensive numerical model for treating surface tension driven flows. As a 
consequence, a full nonlinear model now provides feedback of fundamental interest regarding the appropriate 
interaction of forces for drop formation. The system of numerical programs is lengthy and complex but seems 
adaptable to many problems involving drops and bubbles. 


In a previous study, we considered the so-called "continuous Jet", where equally spaced drops are produced 
in a continuous stream.^ An "intinite Jet" model, assuming periodic boundary conditions for a single drop 
region, provided information on drop stream control through suitable harmonic disturbances. Of general 
interest beyond practical matters of ink Jets is that these numerical solutions gave the complete intervening 
behavior between Rayleigh's Invisc id capillary Jet^ (initial solution) and Lamb's oscillating drop after 
breakups (final solution). 


The present study has the added complication of flow from a nozzle, controlled by a pressure history 
intended for the release of a single drop. A contact circle is now Involved with the uncertainty of what 
constitutes a rigorous treatment of such a boundary. Further uncertainty exists regarding the pressure history, 
produced by a transducer, since the small scale of the Jet precludes measurement of even peak pressures. 
Fortunately the parameter R/W (Reynolds number / Weber number) is in a range that is amenable to numerical 
approximation, and it therefore seems possible that numerical experimentation will shed light on the existing 
uncertainties. 


Currently, the numerical programs are operational with contact circle treatment such that no wetting is 
permitted, as is generally the case with analytic treatment of pendant drops. Motion of the contact circle 
occurs when negative pressures (relative to ambient) are applied at the nozzle entrance and fluid is drawn 
inward from the free surface. Currently, a square wave history of uniform pressure at the nozzle inlet is 
employed to drive the Jet. A positive pulse sufficient to produce the experimental drop size is followed by a 
negative pulse which serves to initiate detachment of a drop. The negative pulse is terminated when the net 
impulse returns to zero. At this stage, surface tension must carry the contraction to final detachment of the 

drop. The objective of the work is to establish quantitatively the conditions which provide optimum drop 

characteristics for high resolution printing. 

In the following, we outline the numerical method which does not differ appreciably from reference 1. Ve 

then discuss a series of results which have thuc far been confined to parameters of interest in ink Jet 

hardware development. Concluding remarks point up those aspects of the work where numerical solution has 
provided practical Insight into the behavior of the flows. 

II. Numerical Method 


The geometry we consider is easily visualized from an examination of output graphics of Fig. 1. A nozzle 
section of unit radius with inlet on the left is assumed to be connected to a cylindrical chamber with radius 
large compared to the nozzle. In Fig. 1 the nozzle section is five units in length. (Note that tick marks on 
plot axes define finite difference mesh distances). Only a small additional annular region is needed in the 
geometry since drops generally will be limited to sizes on the order of the nozzle diameter and wettir.g of the 
outlet facing is not currently permitted. Axi-symmetry is assumed and hence the calculation region involves 
only the upper half of the views in the output graphics. The lower half is added for esthetic purposes in the 
plotting programs. 

The initial solution is impulsive with a prescribed uniform pressure at the circular cross section of the 
inlet. At the initial instant (t«0) the outlet meniscus is flush with the nozzle facing, this being the 
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•qulllbrlum conf Igurition In the ebtence of wetting. Taking the ambient pressure as sero at this free surface 
ve begin with a linear axial pressure gradient in the nossle section. The flow velocities are initially sero. 
The calculation proceeds by explicit tine differences with all succeeding time steps following the same 
prescription which we will now outline for a first step. 


Because we here use a vortlcity streamf unction formulation our first aim is to provide appropriate 
boundary conditions for the latter. This is rather complex since indeed none of the boundaries except the 
axis of symmetry can be treated simply especially since we' want a time varying pressure at the inlet and also 
fish to provide for entrance flow from a chamber, v/e begin in a somewhat roundabout way by computing the 
velocity fields through the momentum equations. In axi-symmetric form these equations are 
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Here u is the axial velocity (z direction) and v ie the radial velocity (r direction). P is the pressure » p 
is Che fluid density and v the fluid viscosity. We choose to express the diffusion terms as vortlcity 
gradients because the vorti- (u)) is always available to simplify the calculation. However « we currently 
are expressing the diff'^sloi* \ the axial velocity equation in terms of u because accuracy in the nozzle 
boundary layer seems to requ^ it. 


A small forward integratiou in tiri. ^f (1) leads to flow because of the pressure gradient in the nozzle. 
Now by integrating the first of (1) over a fluxd (control) volume of the nozzle ve may obtain the uniform 
streamf unction value at the nozzle boundary relative to a zero streamline at the axis. That '<8, the time 
derivative of Lhe streamfunction Q at the nozzle surface becomes 
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vliere we have deleted teims that do not contribute. In (2) L and R refer to left and right cross sections of 

the cylindrical control volume which initially includes the entire fluid region, z is the length of the 

nozzle and r its radius. ° 

o 

We can» of cour.<s, have the control volume continue to be all of the fluid region even when a drop is 
forming but then (2) becomes a more difficult expression. Also, if the free surface draws inward the control 
volume, as expressed in (2), must be shortened by revising z and excluding some of the fluid near the outer 
edge of the outlet. ^ 

Thus far we still do not know the boundary conditions for Q at the entrance of the nozzle or at the free 
surface. For the latter we integrate the expressions 
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singly or in combination, beginning at the axis of symmetry outwx.J to .. i surface. This provides for Q at 
points on the surface at roughly mesh length distances. It is necessary , however, to define the free surface 
as a set of particles generally 5 to 10 per mesh distance. At these surface particles Q is given through 
interpolation among those given above and the value previously obtained at the nozzle boundary. 

Finally, two conditions are provided at the inlet. If flow is Inward ve simulate flow from a large 
chamber by requiring uniform inflow. This effectively requires the boundary layer at the nozzle surface to 
begin at the inlet. With Q already known at the nozzle boundary, a mean velocity u may be obtained using the 
second of (3) for the entire inlet cross section. Reversing the procedure we obtain Q at inlet mesh points. 

Q is not linear in the axis-symmetric case. 

For reverse flow we assume that flow into the chamber will be the same as in the immediately adjacent 
region inside the nozzle or 



there. 

At this point we are not yet ready to solve Internal streamfunction values, we need yet to know the 
vortlcity in the interior. Using the time dependent equation 
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ve IncremAQt th« vorcicity forward In tlno «• we had dona tha valocltiaa. This la of couraa maanlnglaaa for 
the first time stop bacauaa tha vortlclty flald is null. At tha second time step and all following vorticlty 
will derive from the no-sllp nosala surface and to a lesser extant from the free surface. 


Internal streamfunction values are now obtained by simultaneous solution at all net points of the equation 



This is followed by resetting the velocities using (3) and the after thought determination of vorticlty at the 
nozzle surface using 



The redundancy of the above may be questioned but our objective thus far is to provide for more options on 
how solution is to be carried out both for convenience and accuracy. Solution of (5) is of necessity by 
iteration rather than direct because of the time varying fluid region. Convergence is of course enhanced if 
in the course of getting surface streamfunction values through (3) one fills in Internal points. The 
difference in the final outcome from (5) relates to truncation errors in (1), particularly in the nonlinear 
terms. One in fact would find that strict use of (3) only would leave in doubt what value to give the 
streamfunction at the nozzle surface because each integration upward along z"*constant lines would give a 
different error. 

With the streamfunction given at *'tagrangian" free surface points we may evaluate the normal and 
tangential velocities there. That is 
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where T is tangential and n is normal. The second of these is readily obtained from particle Q values • the 
first through a nearest approach method relative to selected interior points followed by interpolation. We 
require surface points or particles to carry with them coordinate information, local axe lengths and angles 
(relative to the axis of sytnnetry). Using the engle we may obtain u and v for the particles through a 
rotational transformation. 

Updating the surface configuration with 


dz 

dt 
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we proceed to determine new local arc lengths along the surface md new angles. Thus 
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It is now easy to obtain surface pressures. Relative to ambient ''zero’* pressure they ere given by 
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«here o is the surface tension coefficient. The first term on the right expresses the destabilizing curvature 
ound the axis of symmetry, the second, the stabilizing curvature In the r-z plane and the last tern the 
pressure contribution from local deformation. For the lest term we again use a closest appr^'^ch aethod for 
evaluation relative to selected Interior points. However, a local continuity expression 
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could probably bo uotd 


At thlo Otago wo may Itorato a Polooon'o aquation for tho prooouro. That io 
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Tho boundary condition for tho prooouro io known ovorywhoro oxcopt at tho notale surface whore the condition 
on tho derivative io 
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G is required at only interior points. 


Obtaining the vorticlty at the free surface has been somewhat of a problem, 
quantities wo write 


In terms of surface 
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Unfortunately in a stationary reference frame this is subject to large errors because in a spherical drop (for 
example) large numbers must effectively cancel to produce a small real vorticlty. To overcome this we define 
velocities relative to a mean for each contiguous part of the fluid so that (14) is evaluated in terms of 
disturbance quantities. 

Scaling to permit nondimensional calculation is the same as with the earlier capillary instability 
problem.^ The reference length is the nozzle radius Tq. The reference velocity is the capillary wave 
velocity 



giving a time scale 



(15) 
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This leads to use of the single parameter, Reynolds number over Weber number 



Thus the above equations a^’e made dimensionless by replacing v by W/R wherever it appears and replacing a by 

1 . 0 . 


III. Results 


During the development of the numerical program there was constant reference to experimental work that was 
being conducted in the development of ink jet hardware.^ Unfortunately because of the small size of the jet 
the crucial measurement of driving pressure could not be made. Lacking this knowledge, along with uncertainty 
concerning the numerical program Itself, made convergence to a successful method very difficult. A failure to 
obtain a solution comparable to experiment could mean being outside the range of pressure values appropriate 
for drop formation or could mean there was a fundamental problem in the numerical program. The program is 
qiiite complex and the running times are by no means short. Unlike experiments which can be conducted in rapid 
succession to achieve a certain operating behavior, the numerical experiments took too much time for many 
wrong guesses. How, however, we arc quite close to experiment and the model has shown evidence of being 
adequate. We will perhaps never know the exact character of the pressure that drives the jet but once having 
achieved agreement in jet characteristics, here demonstrated in a quantitative visual form, we are in a 
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petition to undtrtttnd most proptrtlta of Inctrttt in tht flows. Evtn now, hovtvtr, ont mutt bt mindful of 
llmitttiont of the tpproximtte solution tinet unllmittd grid rttolution it of courts impottibls. 

To illustrtts tht numtrictl solutions vs hsvs chostn s csss which givts tht most inttrtsting ovtrsll 
results schitvtd by the time of this writing. Figures 1» 2 end 3 show s sequence of four solutions etch in 
consecutive times for R/W«5* The pressure history is one which begins with e suction (negative pressure) 
followed by s positive pressure end then e negetive pressure egein. The second negative pressure terminates 
when the net impulse is taro. While not necessary , a very small residual positive pressure rather than sero 
was maintained for the remainder of the running time of this case. The square wave pressure at the inlet and 
times (both in nondimens ional units) are as in Table 1. 


Table 1 


t P/p 


0 - 0.21 -60.0 

0.21 - 0.82 +80.0 

0.82 - 1.43 60.0 

1.43 - + 1.0 


The nondimenslonal times for the plots illustrated are given in Table 2. 

Table 2 


Fig, 

Time 

Plot Interval Q 

QMAX 

QMIN 

Plot Interval P 

PHAX 

PMIN 

1-1 

0.25 

0.03125 

0.0 

-0.52 

4,0 

80,0 

-2.9 

1-2 

0.68 

0.125 

2.03 

0.00 

4.0 

93.3 

0.2 

1-3 

1.03 

0.0625 

1.97 

0.00 

4,0 

7.5 

-84.0 

1-4 

1.64 

0.0625 

0.99 

-0.17 

0,5 

7.1 

-3.0 

2-1 

2.05 

0,03125 

0.78 

-0,08 

0.25 

4.5 

-0.4 

2-2 

2.52 

0.03125 

0.66 

-0.07 

0.25 

4.2 

0.0 

2-3 

2.97 

0.03125 

0.60 

-0.06 

0.125 

4.2 

0.7 

2-4 

3.57 

0.03125 

0.53 

-0.05 

0.25 

7.6 

0.0 

3-1 

4.18 

0.03125 

0.53 

-0.05 

1,0 

17.1 

-1.2 

3-2 

4.52 

0.03125 

0.53 

-0.05 

1.0 

20.7 

-0,4 

3-3 

4.67 

0.03125 

0.54 

-0.05 

2.0 

48.6 

0,6 

3-4 

5.03 

0.03125 

0.52 

-0,03 

2.0 

35.9 

0.6 


In ell figures the streamf unction is plotted as solid lines. Negative streamlines Include tick marks. This 
is manifested as tick marked lines for flow to the left (reverse flow in the nozzle) where the boundary layer 
vorticlty is predominantly negative. The vorticity vector (positive) is outward in upper half plane and 
inward in the lower half. The center line is the streamline Q*0 as discussed previously. Where this line 
appears it is a solid line without tick marks. There are also tick marks on the solid boundaries of the 
nozzle but since these are rectilinear they are distinguishable from streamlines. As mentioned before » the 
tick mark spacing on the rectilinear lines show the grid spacing. This is only roughly true on curved 
streamlines. 

Pressure "isobars" are plotted as dashed lines with shorter dashes for negative pressures. The zero 
pressure line also has long dashes. 

One can use the above information and Table 2 to determine contour values throughout the plots. In some 
places this is difficult but knowing the maximum and minimum values helps and greater detail is probably of 
little use. Because of limited plotting area, say in the contracted neck of the jet, the isobars are 
difficult to distinguish as negative or positive. They are of course positive there because of the high 
surface pressure. Also it should be noted that in the axi-syimetrlc case streamlines are not uniformly spaced 
for a uniform flow velocity. This is the reason for the gap in streamlines along the central axis where a 
tube of flow involves little true flow because its cross sectional area is small, yet the flow velocity may be 
high. 

Beginning with Fig, 1 and referring to Table 1 we note that the time is in the period where the inlet 

pressure is positive but flow is still into the chamber because of fluid inertia from the initial suction 

period. The left most Iwobar at the inlet is +60 (see Table 2) while the lowest pressure is at the concave 

meniscus surface and is -2.9. As a reference it is of interest that, in nondimenslonal units, a drop of unit 

radius has a surface pressure of +2,0. One unit of pressure being being provided by each component of the 
surface tension. Here in the first plot the stable component of the surface tension is such as to restore the 
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flush status of ths nsniscus at tha faes of tha nossla. Acting along with tha posltiva Inlat prassura tha 
flow dlractlon Is soon ravarsad. 

In tha sacond plot of Fig. 1 tha toroidal circulation at tha inlat is baeausa of fl \ banding inward frosi 
tha simulatad larga chambar pracading tha nossla. Tha prassura is nagativa at points just inward fron tha 
circulation whara tha flow rata is high whila at tha cantar lina tha daforaatlon cauaas a prassura axcaading 
tha inlat praraura. Tha maniscus surfaca has a prassura naar unity. 

In tha naxt plot a drop Is forming through forward motion avan though tha inlet prassura has baan 
ravarsad. Nota tha short dashed lines of prassura (nagativa) and a saro prassura in tha contracted region of 
tha incipient drop. Again tha inlat toroidal circulation is still producing a prassura lower than tha applied 
inlat pressure but now more nagativa than tha -60 units there. Tha circulation is in tha process of spreading 
and diaing out. 

In tha last plot of Fig. 1 only remnants of the inlat circulation remains and a separation straamlina now 
divides flow forward in tha drop with predominantly ravarsad flow in tha nossla. Nota that ravarsad flow 
begins in tha boundary layer of tha nossla. At tha time of rha last plot of Fig. 1 tha prassura cycle is 
past. Tha remaining motion of tha Jet is passive except f^r tha minor influence of tha small positive (unit 
prassura) hereafter present at tha inlet. Tha highest prassura is at tha forward most point of tha Jet and is 
a result of local deformations only. 

In Fig. 2 tha progress of tha jet is observed at four succeeding stages. In tha first of these there is 
evidence of local reverse circulation at tha maniscus region tending to remove the concave curvature. This 
and tha continuing flow into tha chamber are because rf tha stabilizing surface tension in this region. At 
the same time the destabilizing surface tension component is causing higher pressures at tha contraction 
although at this stage contraction is dominated simply by mass depletion due to the forward and backward flow. 
These processes continue to the time of the last plot of Fig. 2 whara tha surface tension pressure at tha 
contraction is bi coming an influence on local oahavior. The maniscus region has been pushed outward to be 
convex partly because of the local high pressure and partly from inartial motion from the earlier restoring 
pressure. Also the taper of tha connecting shank to the drop is no longer linear. This further indicates the 
influence of the surface tension pressure at the contraction. 

In Fig. 3 wa continue toward drop break off and Include two plots following breakoff. Note that a 
restoring pressure of the meniscus portion drives fluid into tha chamber to a small extent so that tha outward 
bulge is diminished letvean tha first and second plot of Fig. 3. Also we note that a sacond contraction is 
occurring just behind the drop. 

Drop break off was here preset to occur when only one mesh distance (1/20 of the nozzle diameter) 
remained. At this point the mechanics of the numerical program requires a zero radius and local angles 
consistent with a break. Tha break is assumed to be locally spherical so that the two components of surfaca 
tension arc tha same. This eliminates the singularity of the unstable component which cakes on the value of 
tha stabilizing component. 

Two separate problems are solved simultaneously after breakup. The meniscus converges here to a spherical 
surface of unit radius because of tha residual unit prassura at tha inlet. The drop portion of the Jet 
proceeds toward a state where a drop plus a satellite exist because of tha contraction behind tha main drop. 

IV. Concluding Comments 

Wa are interested here in summarising those things learned from tha numerical model thus far. In soma 
cases the Information was partially known from experiment and matters ware clarified and/or confirmed by tha 
model. In other cases tha information from the model gave insight thus far unobtainable experimentally. 

1. The peak pressure magnitudes for driving tha Jet are now obtainable. These » of course, depend on nozzle 
losses but can be fairly well estimated for an initial trial calculation. 

2. The pressure history to drive the jet is nearly a symmetric one for clean drop formation with a net saro 
impulse for a single drop. The negative phase of the prassura history is essential to drop break off in 
the range of parameters here studied. 

3. Uniform prassura at the inlat cross section and ignored watting of tha outer facing of the nozzle are 
reasonable simplifying assumptions that do not affect comparison with experiment. 

4. In tha range of parameters hare studied, the meniscus is vary ovardampad as is tha drop. Experimental 
maniscus oscillations presumably relate to chamber or chambar to nozzle transition characteristics. Tha 
latter can be studied with tha present modal. 

5. Tha flow fields throughout tha history of drop generation may be analyzed in almost complete detail from 
tha numerical results and, internal flows, usually not visible experimentally, may hare be studied. 
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